Ostwald ripening of faceted two-dimensional islands 
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We study Ostwald ripening of two-dimensional adatom and advacancy islands on a crystal surface 
by means of kinetic Monte Carlo simulations. At large bond energies the islands are square-shaped, 
which qualitatively changes the coarsening kinetics. The Gibbs-Thomson chemical potential is 
violated: the coarsening proceeds through a sequence of 'magic' sizes corresponding to square or 
rectangular islands. The coarsening becomes attachment-limited, but Wagner's asymptotic law 
is reached after a very long transient time. The unusual coarsening kinetics obtained in Monte 
Carlo simulations are well described by the Becker-Doring equations of nucleation kinetics. These 
equations can be applied to a wide range of coarsening problems. 
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I. INTRODUCTION 

Domains of a guest phase inside a matrix tend to 
coarsen, thus reducing their specific interface energy. The 
prominent mechanism of coarsening was proposed by 
Ostwald 1 more than hundred years ago: larger domains 
grow at the expense of smaller ones by exchanging atoms. 
The net atom flux is directed to larger domains since they 
possess smaller interface energy per atom. The semi- 
nal theory of Ostwald ripening was proposed by Lifshitz 
and Slyozov^ and by Wagner 4 They showed that, at late 
times, the system is characterized by a single characteris- 
tic scale, namely, the average domain size R(t). The time 
evolution of the system consists in changing the scale: the 
domain distribution, shape of the diffraction peaks, etc. 
remain unchanged when scaled by R(t). The average do- 
main size follows, in turn, universal laws, R(t) cx i 1 / 3 
if the atom diffusion is the rate limiting process^ and 
R(t) oc t 1 ' 2 if the attachment-detachment at the domain 
interface is the limiting onei 3 - 

The kinetic scaling is essentially based on the Gibbs- 
Thomson formula /i = 7 / R for the excess chemical po- 
tential of a gas that is in equilibrium at the curved sur- 
face of a liquid droplet (the constant 7 is proportional 
to the surface tension). The aim of the present work 
is to study the Ostwald ripening kinetics at low tem- 
peratures (or large bond energies) when the crystalline 
droplets are faceted. The energy of a small crystalline 
droplet is minimum at 'magic' sizes when all facets are 
completed. The coarsening proceeds as a sequence of 
jumps from one magic size to the next. We perform ki- 
netic Monte Carlo simulations of Ostwald ripening ki- 
netics for faceted two-dimensional (2D) islands and find 
very long transient behavior of the system, so that the 
universal asymptotic laws are still not reached. We de- 
velop a mean-field theory for Ostwald ripening, based on 
the Becker-Doring 4 equations. We show that these equa- 



tions, being the basic equations of nucleation theory;^ 
can be used to describe the coarsening kinetics in the 
whole size range, starting from monomers up to the long- 
time asymptotics that are not available in Monte Carlo 
simulations. Both the Lifshitz-Slyozov- Wagner regime 
and the coarsening through a sequence of magic sizes are 
well described. This approach requires only the knowl- 
edge of the droplet energy dependence on the number of 
atoms in the droplet and can be applied to a wide range 
of coarsening problems in other systems as well. 

Two-dimensional (2D) islands on a crystal surface are a 
practically important physical system that reveals differ- 
ent coarsening mechanisms and allows detailed theoret- 
ical and experimental studies of the coarsening kinetics. 
From the experimental studies, we mention the ones that 
report time exponents n in the coarsening law R(t) cx 
t n . These include low-energy electron diffraction from a 
chemisorbed monolayer of oxygen on W(110),— ^ helium 
atom beam diffraction from 0.5 monolayer (ML) of Cu on 
Cu(100), 9 optical microscopy of a thin layer of succinon- 
itrile within the liquid-solid coexistence regio n 10 ' 11 and 
a binary mixture of amphiphilic molecules, 1 ?, and low- 
energy electron microscopy of Si on Si(001) i 13 i 14 In these 
works ) 7 ' 8 ' 9 ' 10 ' 11 ' 12 the time exponents somewhat smaller 
than 1/3 were found and explained by the Lifshitz 
Slyozov law with finite-size corrections. The time ex- 
ponent 1/2 obtained for Si on Si(001)i2J^ was treated as 
the case of kinetics limited by the attachment and detach- 
ment of adatoms to steps£ Our recent x-ray diffraction 
study of coarsening of 2D GaAs islands on GaAs(OOl), 15 
which showed an apparent time exponent close to 1, was 
the experimental inspiration for the present work. 

Two-dimensional islands of 'magic' sizes were observed 
on several surfaces, such as Pt(lll)^, Si(lll)r 7 - and 
Ag(lll)^ (see also a review^). It was shown theoret- 
ically that the presence of magic island sizes disrupts 
the scaling law of submonolayer molecular beam epitaxy 
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growth. 20 Magic sizes of three-dimensional Pb nanocrys- 
tals on Si(lll) lead to a breakdown of the classical Ost- 
wald ripening laws. 21 

Monte Carlo simulations of Ostwald ripening were per- 
formed using the 2D Ising mode L 22 ' 23 ' 24 They were lim- 
ited to rather small values of the coupling constant, so 
that the domains are rounded and faceting is absent. The 
time exponents were found to be smaller than 1/3, which 
was explained by finite-size corrections to the Lifshitz- 
Slyozov law. Further discussion of theoretical and simu- 
lation studies can be found in several reviews ! 25 ' 26 ' 27 De- 
spite kinetic Monte Carlo simulations are routinely used 
to model epitaxial growth ) 28 ' 29 ' 30 ' 31 ' 32 ' 33 we are aware of 
only one such study of coarsening of 2D islands on a 
crystal surface.— This latter simulation was limited to 
small bond energies and rounded islands, similar to the 
simulations of the Ising model. 

A physical difference between the coarsening kinetics 
of 2D epitaxial islands and that of Ising spins becomes ev- 
ident when we compare adatoms and advacancies on one 
side with up and down spins on the other side. The first 
two objects possess qualitatively different kinetics (mo- 
tion of an advacancy is a result of the collective motion 
of atoms), while up and down spins are equivalent. This 
distinction manifests itself in the transition probabilities, 
as discussed below. The fundamental laws of Ostwald 
ripening are expected to be independent of the tran- 
sition probability distribution, so that a kinetic Monte 
Carlo simulation of the coarsening of epitaxial islands 
allows one to check this conclusion. Here, we perform ki- 
netic Monte Carlo simulations of Ostwald ripening of 2D 
adatom islands (surface coverage 0.1 ML) and 2D adva- 
cancy islands (surface coverage 0.9 ML) in a wide range 
of bond energies (or temperatures) . Our particular aim is 
to perform simulations in the case of large bond energies 
(low temperatures) when the islands are faceted, which 
was not studied previously. 



II. MONTE CARLO SIMULATIONS 
A. Simulation method 

We employ the well-established generic model devel- 
oped for kinetic Monte Carlo simulations of molecular 
beam epitaxyj 28 i 29 i 30 i 31 i 32 i 33 i 34 Atoms occupy a simple cu- 
bic lattice and interact with a pair energy that depends 
only on the number of bonds. An alternative approach 
to simulate surface kinetics is a detailed Monte Carlo 
simulation of a particular surface with energetic param- 
eters taken from ab initio calculations, as it was done 
for GaAs(OOl) or InAs(001)^i 3 ^i 3 ii2 8 .^ Such simulations 
are very time-consuming and hence are limited to small 
time and spatial scales. They can hardly be applied to 
study the coarsening process. Some characteristic fea- 
tures of compound semiconductors can, however, be in- 
cluded in the generic model as a compromise! 40 ! 41 ' 42 

We use an algorithm 4 — that advances simulated time 



depending on the probability of the chosen event. This 
algorithm is commonly used in the epitaxial growth simu- 
lations. We note that the Ostwald ripening simulations of 
the 2D Ising mode l 22 ' 23 ' 24 have employed the Metropolis 
accept reject algorithm. This algorithm becomes ineffec- 
tive at low temperatures, since most of the attempts are 
rejected and computer time is wasted. That is why pre- 
vious simulation s 22 ' 23 ' 24 were restricted to relatively high 
temperatures T > 0.5T C , where T c is the Ising phase tran- 
sition temperature. Of course, both algorithms give the 
same results and differ only in the computation time. 

The choice of the probability w(x — > y) for the tran- 
sition from the state x to the state y incorporates the 
physics of the system into the simulations. The choice 
is made differently for the epitaxial growth and the 
Ising model simulations. It is worthwhile to compare 
these probabilites briefly. A sufficient condition that 
the system evolves to thermodynamic equilibrium is the 
detailed balance condition, w(x — > y)/w(y — > x) = 
exp(-AE/k B T). Here AE = E(y) - E(x) is the energy 
difference between the states x and y, &b is the Boltz- 
mann constant and T is the temperature. The simula- 
tions of the Ising model use a probability that depends on 
AE (either the Metropolis or the Glauber probability). 
These probabilities favor transitions which reduce the en- 
ergy of the system, AE < 0. On the other hand, for an 
atom jump on the crystal surface, the transition proba- 
bility does not depend on the final state y but only on the 
height of the energy barrier that needs to be overcome^ 
The probability is w(x — » y) oc exp[E(x)/kBT], where 
E(x) < is the energy of the initial state with respect 
to the barrier. Such a probability obviously satisfies the 
detailed balance condition. The system evolves into a 
lower-energy state since it escapes higher-energy initial 
states with larger probabilities. 

In the present study, no step edge barrier is imposed. 
An atom detaching from a step edge can go to the lower 
or the upper terrace with equal probabilities. In particu- 
lar, atom exchange between advacancy islands is achieved 
predominantly by adatoms diffusing on the top level 
rather than by the diffusion of vacancies, despite that the 
latter process is not forbidden. Similar simulations, but 
with an infinite step edge barrier, were performed in our 
preceding work4£ In that study, the restriction for atoms 
to escape a vacancy island to the higher level resulted in 
another coarsening mechanism, diffusion and coalescence 
of whole islands due to atom detachment and reattach- 
ment within an island. The coarsening by dynamic coa- 
lescence is much less effective than Ostwald ripening and 
becomes essential when the detachment of atoms form 
islands is prohibited. 

An atom that has n neighbors in the initial state 
with equal bond energies Eb to these neighbors, pos- 
sesses an energy E(x) = —(nEt, + Ed), where the ac- 
tivation energy of surface diffusion Ed is the barrier 
height. It determines the time scale r of the problem, 
t -1 = j/exp(— En/k^T), where v w 10 13 s _1 is the vi- 
brational frequency of atoms in a crystal. In the epitaxial 
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growth simulations, the time scale r is to be compared 
with the deposition flux, which determines an appropri- 
ate choice of Ed- We do not consider deposition, and 
the choice of Ed is arbitrary. Note that the works on 
the Ising model kinetics measure time simply in the flip 
attempts (sweeps) per lattice site. We take the same val- 
ues of Ed as in the preceding work^ with the aim to 
compare time scales of Ostwald ripening (in absence of 
the step edge barrier) with that of dynamic coalescence 
(infinite step edge barrier). Namely, we take Ed = 0.2; 
0.1; eV for E b = 0.2; 0.3; 0.4 eV, respectively. 

The ratio of the interaction energy between neighbor- 
ing atoms to the temperature E^/k^T is the only es- 
sential parameter for the coarsening problem. We fix 
the temperature at 400 K and vary the bond energy Eb 
from 0.2 eV to 0.4 eV. In terms of our model, the Ising 
phase transition takes place at E^jk^T = 21n(l + \/2). 
Our choice of bond energies corresponds to T/T c varying 
from 0.15 to 0.3, temperatures much lower than the ones 
used in previous kinetic Monte Carlo studies of Ostwald 
ripeningi 22 i 23 ' 24 i 34 Here T c is the Ising phase transition 
temperature. 

We perform kinetic Monte Carlo simulations on a 
1000x1000 square grid with periodic boundary condi- 
tions. Each simulation is repeated 25 times, to obtain 
sufficient statistics for the island size distribution. In 
the initial state, either 0.1 ML or 0.9 ML are randomly 
deposited. Adatom islands form in the first case and ad- 
vacancy islands in the second. 



B. Simulation results 

Snapshots of the simulated system at the end of a sim- 
ulation are presented in Fig.QJa). As the bond energy Ef, 
is increased (from left to right), the island shape contin- 
uously transforms from more circular to almost square. 
Since faceting transitions are absent in 2D systems, we 
refer to the almost square islands as faceted in order to 
stress the qualitative shape difference at small and large 
bond energies. Apart from the change in shape, the equi- 
librium density of adatoms between islands exponentially 
decreases as the bond energy increases. 

Figures [1Kb) an d (c) show time variations of average 
island diameters 2R(t) in logarithmic and linear scales, 
respectively. The determination of an average island size 
is described in Sec. Ill Cl At small bond energies (left col- 
umn in Fig. [1]), the process of Ostwald ripening follows 
the Lifshitz-Slyozov law R(t) cx t 1 ^ 3 . As the bond energy 
increases, the coarsening law for advacancy islands devi- 
ates from that for adatom islands and from the expected 
t 1 / 3 law. At large bond energies (right column in Fig.[T]), 
the coarsening behavior of advacancy islands is qualita- 
tively different and close to a linear dependence, in a wide 
range of island sizes. The coarsening of adatom islands 
also notably deviates from the Lifshitz-Slyozov law. The 
attachment-limited asymptotic t 1 ^ 2 can be inferred from 
the figure, but it is not really reached. 



Figure QJd) shows the island size distributions at dif- 
ferent times. The uniformly spaced time instances are 
marked on the curves in Fig.rjjc) by the same symbols as 
used for the corresponding size distributions. The distri- 
butions are scaled by the average size R(t) : instead of the 
probability P(r), we plot RP(r) versus r/R. The scaled 
distributions do not change in time even at large bond 
energies, where the average island sizes do not show a 
power law behavior. The island size distribution notably 
changes with increasing bond energy, Fig.[T]Jd). The dis- 
tribution develops a tail extended to 2i?, while at smaller 
bond energies it is limited to 1.5R. 



C. Analysis of the simulation data 

We obtain the sizes of all islands in the simulated sys- 
tem by using an algorithmic that allows to count all topo- 
logically connected clusters in the system. At large bond 
energies, we average the radii r n = yripn (where n is 
the number of atoms in a cluster) of all islands, exclud- 
ing individual adatoms from the distribution. In the case 
of small bond energies we find that, besides monomers, 
a notable amounts of transient dimers, trimers, etc. are 
present in the simulated system. Their densities quickly 
decrease with increasing number of atoms in the clus- 
ter and they are well separated from the distribution of 
the large clusters. If these small clusters are included in 
the island size distribution when calculating average ra- 
dius R, we obtain unreasonable time dependencies R(t). 
Hence, we calculate the averages taking into account is- 
lands of at least 6 atoms. 

We also use the Monte Carlo simulations to verify the 
average island size determination in diffraction studies. 
In a diffraction experiment, one has access to the peak 
profile only and obtains the average size from its width. 
Using the island distribution obtained in the simulation 
and calculating the peak profiles, we can compare the 
average sizes obtained from the real space and the recip- 
rocal space distributions. The diffraction peaks (struc- 
ture factors) obtained from the simulations are present 
in Fig. (Ha). We consider the anti-Bragg condition (sub- 
sequent atomic layers contribute to the scattering func- 
tion with a phase shift of 7r) and obtain two-dimensional 
intensity distributions I(q x ,q y ) from Fourier transforma- 
tion of exp[iirh(x,y)]. Here an integer function h(x,y) is 
the surface height. Then, we take into account that in 
a diffraction experiment, the scattered intensity is usu- 
ally collected by a wide open detector that integrates 
over one of the components of the scattering vector 
Hence, we integrate the distributions I(q x ,Qy) over one 
of the components of the scattering vector, either q x or 
q y . The resulting diffraction peaks I(q) are presented in 
Fig. HJa). The peaks corresponding to different time mo- 
ments [the same time moments as in Fig. [ljd)] coincide 
after the wave vectors q are scaled by the average island 
size. Kinetic scaling is thus confirmed. The shapes of the 
peaks depend on the bond energy E^, thus showing that 
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FIG. 1: Results of kinetic Monte Carlo simulations: (a) snapshots of the 1000X 1000 simulation cells at the end of the simulations, 
(b) and (c) time dependence of the average island size in logarithmic and linear scales, and (d) the island size distributions. 
The gray levels in the snapshots vary from black to white as the surface height increases. Different columns show results for 
different bond energies Eb, with the temperature fixed at T = 400 K. The size distributions are obtained at the time moments 
marked in (c) by the corresponding symbols. 



the island size distribution and the correlations between islands change. 
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FIG. 2: (a) Diffraction peaks obtained from the Monte Carlo 
simulation results (the gray curves are Gaussian fits), (b) 
Time dependencies of the average island sizes obtained from 
the numbers of atoms in the clusters (black curves) and from 
the widths of the diffraction peaks (gray curves). 



The quantity most commonly measured in a diffraction 
experiment is the full width at half maximum (FWHM) 
of a peak obtained by an appropriate fit. Considering 
islands of linear size 2R, one obtains a structure fac- 
tor sin 2 (gi?)/sin 2 (qa), which can be approximated by 
exp(— q 2 R 2 /n)^ Here, a is the lattice spacing. We ob- 
tain the average size 2R by fitting the peaks to this Gaus- 
sian function, despite the peaks are not Gaussian, es- 
pecially for small bond energies. Figure [Hb) compares 
these sizes with the ones obtained from the real-space 
island size analysis described above. The values are in 
good agreement, thus confirming that the average quan- 
tities can be obtained from the diffraction peak widths 
even if the profiles deviate notably from Gaussian. 



III. COARSENING EQUATIONS 

A. The Becker— Doring equations for the 3D 
problem 

The process of Ostwald ripening can be described by 
two alternative approaches, either in terms of a contin- 
uous function f(r) representing the number density of 
clusters of radius r, or in terms of discrete numbers c n 
representing the densities of clusters containing n atoms 
(nmers). The first approach was employed by Lifshitz 
and Slyozov^ and Wagner^ The equations for discrete 
quantities c„ were first formulated by Becker and Doring 4 
and ever since form the basis of nucleation theory^* 6 - 
Closely related equations, the rate equations, were used 
in the description of crystal growthj 48 i 49 i 50 They contain 
an additional deposition term, while the detachment pro- 



cess is not essential and the corresponding terms in the 
equations are frequently omitted. Similar discrete equa- 
tions for the Ostwald ripening process were introduced 
under the names of microscopic continuity equations ; 51 ! 52 
population balance equations ; 53 ! 54 ' 55 or rate equation 
approach. 56 Mathematical aspects of the relationship be- 
tween the discrete and the continuous equations were also 
considered j 57 ' 58 The aim of the present section is to link 
the discrete and continuous approaches and obtain equa- 
tions that can be used for a numerical study of the Ost- 
wald ripening process. 

The number of atoms n in a cluster increases or de- 
creases by one when an atom is attached to the cluster 
or detached from it. Let J„ be the net rate of transfor- 
mation of nmers into (n + l)mers. The number c„ of 
nmers increases due to the transformation of (n — l)mers 
into nmers and decreases because of the transformation 
of nmers into (n + l)mers: 

dc n /dt = J n _i - J„. (1) 

This equation is valid for n > 2. The equation describing 
the number of monomers c\ is obtained by requiring that 
the total number of atoms in the system 

oo 

N=J2 nC n (2) 

n=l 

does not change in time. The condition dN/dt = gives, 
after substitution of Eqs. |T]) and rearrangement of the 
terms, 

oo 

da/dt = -2Ji - Jn- (3) 

n=2 

This equation takes into account that each transforma- 
tion of an nmer into an (n + l)mer decreases the number 
of monomers by one, except in the case n — 1, where two 
monomers form a dimer. 

The net rate J n is a result of two processes. First, 
an nmer catches a monomer. The rate of this process 
is proportional to the densities of the nmers and the 
monomers and can be written as & n c\c ni where a n is 
a time-independent coefficient that remains to be deter- 
mined. The second process is a spontaneous detachment 
of a monomer from a (n + l)mer. It is proportional to 
the density of (n + l)mers solely and can be written as 
b n c n +ii where b n is another time-independent coefficient 
to be specified. Hence, we obtain 

J n = a n c x c n - b n c n+1 . (4) 

Equations (P), ([3]), and ^ are the Becker-Doring equa- 
tions. 

If the time limiting process is the adatom diffusion be- 
tween clusters, the attachment and detachment coeffi- 
cients a n and b n for the 3D problem are calculated, for 
large n, as follows. The cluster of n atoms is consid- 
ered as a sphere of radius r n , so that n — Airr^/S. To 
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calculate the attachment coefficient, we solve the steady- 
state diffusion equation V 2 c(r) = with two boundary 
conditions: the concentration of the monomers far away 
from the cluster is equal to their mean concentration, 
c(r) \ r =oo — ci, while the concentration of the monomers 
at the cluster surface is zero, c(r) |r=r n = 0, since the 
monomers are attached to the cluster as soon as they 
reach it. The solution is c(r) = (1 — r n jr)c\. The total 
atom flux at the cluster surface 

j n = 4:TTrlDVc(r) | r=r „ , (5) 

where D is the diffusion coefficient of the monomers, is 
equal to A / nDr n c\ 1 and hence the attachment coefficient 
is 

a„ = AitDr n . (6) 

The detachment coefficient is calculated assuming that 
the concentration of the monomers at the cluster sur- 
face is equal to the equilibrium monomer concentra- 
tion c ncq , while there is an ideal sink for monomers at 
infinity, c(r) \ r=oc = 0. The solution of the steady- 
state diffusion equation with these boundary conditions 
is c(r) = c ncq r„/r, and the corresponding detachment 
flux of the monomers is 6„+i = 47rDr n c neq . Here we 
take into account that this flux refers to the detachment 
from the (n + l)mer. The ratio of the detachment and 
the attachment coefficients is then 

(7) 

The equilibrium density of monomers at the surface of a 
cluster is given by the Gibbs-Thomson formula 

c„ C q = c ooeq exp(7/r„) w c oocq (l +j/r n ), (8) 

where 7 is a constant proportional to the surface tension. 
The explicit expression for 7 is given in the next section. 
A correction to Eq. {8| for small clusters consisting of 
very few atoms, that is important in the theory of nuclc- 
ation, is not essential for the Ostwald ripening problem. 
Then, equations CD)-© give a complete set of equations 
that describe the process of Ostwald ripening. 

When clusters are large enough, n can be treated as 
a continuous variable. Let us verify that the continuous 
equations derived from the set of equations above are the 
Lifshitz-Slyozov equations. The cluster size distribution 
function f(r,t) is defined so that f(r,t)dr is the num- 
ber of clusters per unit volume in an interval from r to 
r + dr. Then, f(r,t)dr — c n (t)dn and, keeping in mind 
that n = 47rr 3 /3, we obtain f(r, t) = Anr 2 c n (t). The 
mass conservation law © can be rewritten, by separat- 
ing monomers and larger clusters, as 

47T f°° 

c iW + ^-/ r 3 f(r,t)dr = N = const. (9) 
3 Jo 

The finite-difference equation |T]) transforms into the con- 
tinuity equation 

df/dt + dJ/dr = Q. (10) 



To calculate the flux in the cluster size space J(r,t), one 
can neglect the difference between c„ and c n+ i in Eq. (f?|). 
Then, substituting Eqs. ([7]) and {H}, one obtains 

j(r i t) = -(c 1 -c ooe(l -^l)f. (11) 



Equations (f9"|- (|ll[) coincide with the Lifshitz-Slyozov 
equations £ 




FIG. 3: (a) The cluster size distribution obtained by numeri- 
cal solution of the Becker-Doring equations at different times 
(thin black lines, the lines closer to the gray line correspond to 
later times) and the analytical solution by Lifshitz and Sly- 
ozov (thick gray line), (b) The time dependency of R 3 . A 
linear asymptotic is evident from the plot. 

As an example, we compare in Fig. [3] numerical solu- 
tions of the ordinary differential equations (HJ-© with 
the analytical result^ To solve the Becker-Doring sys- 
tem, we employ a second-order Rosenbrock method, 
which is essentially based on a Pade-approximation of 
the transition operator (see, e.g., Ref. [591 ). A version of 
this method^ that fits well to stiff systems of differential- 
algebraic equations was used. Practically, we solve a set 
of up to one million ordinary differential equations on 
a personal computer. The solutions in Fig. [3] are ob- 
tained by taking 7 = 5 and, as the initial condition at 
t = 0, only monomers with the initial supersaturation 
ci/coocq = 10 5 . The figure shows that the numerical so- 
lutions asymptotically converge to the analytical formula, 
which validates our approach. 

B. Attachment and detachment coefficients 

Equation |7]) can be derived in a more general form 
that will be useful for the considerations below. In equi- 
librium, all fluxes J n are identically equal to zero. Then, 
denoting by C n the equilibrium concentrations of the 
nmers, we have from Eq. ((4|) 

b n +i/a n = CiC n /C n +i. (12) 

The equilibrium concentrations calculated in the frame- 
work of equilibrium thermodynamics are 61 

C n = CI exp[-(£ n - nE^/k^T], (13) 

where E n is the energy of an nmer and E\ is the en- 
ergy of a monomer. This relation can be treated as the 
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mass action law for the equilibrium between nmers and 
monomers, C n ^ nC\. Substitution into Eq. (|12p gives 



b n+ i/a n = Coo Cq exp[(i;„ + i - E n )/k B T] 



(14) 



where Coo C q = exp(— E\/kT) is the concentration of 
monomers that are in equilibrium with an infinite cluster. 
For spherical clusters, Eq. (fT4|) reduces to the Gibbs- 
Thomson formula. The energy of a spherical cluster is 
E n = 4irr 2 a, where a is the surface tension, with the 
radius r defined by nv = 4-7rr 3 /3, where v — a 3 is the 
volume per atom. The radius increase due to the attach- 
ment of an atom to a nmer is given by v — 4irr 2 Ar. The 
change of the energy due to the attachment of a single 
atom is E n+ \ — E n — 8tt or Ar = 2va/r. Thus, we arrive 
at Eq. (JSj) with 7 = 2vo/k-oT. A similar calculation for 
the 2D case gives 7 = so /k^T, where s is the area per 
atom. 
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FIG. 4: The island size distribution of faceted islands ob- 
tained in the kinetic Monte Carlo simulations (a) and by nu- 
merical solution of the Becker-Doring equations (b). The 
strong preference of magic island sizes is obvious. 

Equation (fl"4"l) is more general than the Gibbs- 
Thomson formula and can be used in situations when 
the latter is not applicable. Figure SJa) presents the is- 
land size distribution obtained in our kinetic Monte Carlo 
simulations at an early stage of coarsening for the largest 
bond energy we have studied, Eb = 0.4 eV. The distribu- 
tion is not smooth but consists of peaks at 'magic' island 
sizes corresponding to a product of two close integers, 
like 30 = 6 x 5. Accordingly, the insert in the figure 



shows that the islands are mainly rectangles with an as- 
pect ratio close to 1. The origin of such a distribution 
is evident: when an island consisting, for example, of 30 
atoms, grows by one atom, its energy increases by 2Eb, 
while further growth to 36 atoms does not change its en- 
ergy at all. Thus, we solve the Becker-Doring equations 
with the energy of a 2D island of n atoms calculated 
as follows. First, we find the largest square that still 
contains fewer atoms than n. Then, we add, as long as 
the number of atoms does not exceed n, rows of atoms 
to the side of the square. The last row may be incom- 
plete. The number of broken bonds for such an island is 
calculated. Figure HJb) presents a numerical solution of 
the Becker-Doring equations with the island energies E n 
thus calculated and the attachment-detachment coeffi- 
cient ratio given by Eq. (|14[) . The approximation for a n 
appropriate for the 2D case is given below in Sec. IIII CI 
The size distribution closely reproduces the one obtained 
in the Monte Carlo simulations: squared or rectangular 
(with aspect ratio close to 1) islands are discrete barriers 
to be overcome, while the filling of an atomic row does 
not change the island energy and proceeds relatively fast. 
This example shows that Eq. (fT4"]) can be used when the 
island energy E n is known but is not described simply by 
the surface tension, and the Gibbs-Thomson formula is 
not applicable. 



C. Coarsening equations in two dimensions 

The Becker-Doring equations (HJ)— (J3J) and the equa- 
tion (fT4f for the ratio of the coefficients b n +i/a n do not 
depend on the dimensionality of the system and can be 
applied to both 2D and 3D problems. (It may be worth 
to note that the radius r n entering the Gibbs-Thomson 
law is expressed differently through n in the 2D and 
3D cases.) The only formula that has to be reconsid- 
ered is expression ([5]) for the attachment coefficients a n , 
since it is based on the solution of the 3D diffusion equa- 
tion. The solution of the 2D diffusion equation behaves 
as c(r) cx lnr and the boundary condition c(r) | r=00 = c\ 
cannot be imposed. A simple approximation is to place 
this condition at a finite distance Z, given by an aver- 
age distance between the islandsi 52 i 62 ' 63 i 64 ' 65 Then, in the 
case of diffusion-limited kinetics, the attachment coeffi- 
cient a n does not depend on n and is proportional to 
(InZ) -1 . Proceeding to the continuous distribution func- 
tion, one arrives at Eq. (jlip . with the conservation law 
© rewritten for the 2D case. The coarsening equations 
are solved analytically in this case i 63 i 64 ' 66 

A self-consistent description of two-dimensional diffu- 
sion can be obtained by taking into account its screen- 
ing by the island distribution*^ A solution of the 2D 
screened diffusion equation, satisfying the boundary con- 
ditions c(r) | r=00 = ci and c(r) \ r=r „ = 0, is c(r) = 
ci[\ — if (r/£)/-Ko( r ri/OL where K (x) is the zeroth 
modified Bessel function and £ is the screening length 
that remains to be defined. Then, one obtains the at- 
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tachment coefficient 

a n = D/C(r n /0, (15) 

where 

K(x) = 2irxK 1 (x)/K (x) (16) 

and K\(x) is the first modified Bessel function. The self- 
consistency condition for the screening length £ is 67 

/>oo 

r l = / lC(r/Of(r,t)dr. (17) 
Jo 

Expressions very similar to Eqs. (f!5]) and (fl6|) are used 
in studies of crystal growth from the gas phase d 49 ' 50 , 
with one essential difference: for the latter problem, the 
length £ is the mean diffusion length of an adatom on 
the surface before its reevaporation. It is a well-defined 
time- independent constant, so that no self-consistency 
condition is involved. 

In the case of attachment-limited kinetics, the bound- 
ary condition for the concentration field c(r) at the is- 
land surface is the absence of the flux, Vc| r=r „ = 0, 
which gives a constant solution, c(r) = c%. Then, the 
attachment coefficient is 

a„ = 2nKr n , (18) 

where K is the attachment coefficient. The result is in- 
dependent of screening effect in this case. The same ex- 
pression is obtained in the approximation of a constant 
screening distance equal to the mean distance between 
islands, 52 i 62 i 63 i 64 i 65 



D. Coarsening equations for advacancy islands 

In our Monte Carlo simulations, a step edge barrier 
is absent and an atom detaching from a vacancy island 
ascends to the higher terrace. The vacancy island size 
increases by a vacancy at the same time. The coarsen- 
ing proceeds by exchange of adatoms between vacancy 
islands and can be described by equations similar to the 
Becker-Doring equations. Let us denote by g(t) the con- 
centration of adatoms, while c n are the concentrations of 
2D islands of n vacancies. Then, the continuity equa- 
tion (fTJ) for the density of clusters c n (t) remains un- 
changed. The fluxes J n in these equations describe two 
processes. The first process is the spontaneous emission 
of an adatom. Its rate is proportional to the density 
of nmers. The second process is an absorption of an 
adatom by the vacancy type (n + l)mer, which gives rise 
to a fimer. Its rate is proportional to the density g of 
adatoms and the density of (n + l)mers, so that 

J n = b n c n - a n+ igc n+ i. (19) 

The annihilation of an atom and a single vacancy is 
described by the flux Jq — —a\gc\. Then, the set of 



equations (Q} is valid for n > 1. The creation of an 
adatom-vacancy pair from a flat surface is prohibited in 
our model. 

Since the growth of a vacancy cluster by one vacancy is 
accompanied with the emission of one adatom, the con- 
served total amount of atoms in the system is given by 

oo 

N = ^ nJ n — g, (20) 

n=l 

which replaces Eq. @. By differentiating this equation 
with respect to time and rearranging the terms, we obtain 
from dN/dt = an equation for the time variation of the 
adatom density: 

oo 

dg/dt = J2 Jn- (21) 

n=0 

The mass action law now has to be written for an equi- 
librium between an advacancy island and adatoms that 
annihilate, C n + ng ^ 0. Hence, instead of Eq. (fT"3f we 
have 

C n g n = exp[-(£„ + nE x )/k B T}. (22) 

The requirement of zero fluxes at equilibrium gives rise 
to the detailed balance condition 

b n /a n+1 = Coo 0q exp[- (E n+1 - E n )/k B T] (23) 

that differs from Eq. (TTJJ by the sign in the exponent. For 
circular islands, the same calculation as above leads to 
the Gibbs-Thomson formula ([5]) with negative 7, which 
corresponds to a negative curvature of the vacancy island 
surface. 



E. Solutions of the coarsening equations 

Figure [5] presents the results of the numerical solution 
of the Becker-Doring equations for adatom and adva- 
cancy islands. With the aim to quantitatively compare 
the solutions with the results of kinetic Monte Carlo sim- 
ulations in the whole time interval, we obtain the average 
island sizes in the same way as in the simulations, and use 
the same initial conditions. Namely, the average island 
sizes are calculated taking into account the islands con- 
taining at least 6 atoms, for the reasons discussed in Sec. 
Ill CI The initial random adatom distribution with the 
coverage 0.1 ML contains not only monomers, but also 
dimers, trimers, etc., the densities of which quickly de- 
crease with the increasing number of atoms in the cluster. 
By simple statistical analysis of the initial distribution in 
kinetic Monte Carlo simulations, we find that at t = 0, 
c„ ps ci x lO'" -1 )/ 2 . This distribution was used as the 
initial condition for the Becker-Doring equations. The 
initial conditions are essential only at the initial stages 
of coarsening. The results of the calculations do not de- 
pend on the initial monomer concentration c\ , as long as 
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time (s) time (s) 
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r/R r/R 

FIG. 5: The results of numerical solution of the Becker- 
Doring equations: time dependencies of the average island 
sizes in logarithmic (a) and linear (b) scales and the island 
size distributions (c). The left column presents calculations 
for the bond energy Eb = 0.2 eV with diffusion-limited ki- 
netics, while the right column shows the results for the bond 
energy Eb = 0.42 eV with attachment-limited kinetics. The 
solutions of the Becker-Doring equations are shown in (a) and 
(b) by black lines, and the results of the kinetic Monte Carlo 
simulations by gray lines. Symbols "a" and "v" on the plots 
denote the results for adatom and advacancy islands, respec- 
tively. Full lines in the right column show the calculations for 
the discrete island energies with 'magic' sizes, while the bro- 
ken lines are calculations for the continuous Gibbs- Thomson 
chemical potential. 



the initial supersaturation c\ (t — 0) /coo C q is much larger 
than unity. The time scales of the solutions are adjusted 
to these of the Monte Carlo simulations. 

The case of small bond energies (left column in Fig. 
[5]) is well described by the 2D diffusion limited kinetics 
with screening (| 1 5|) and the ratio of the detachment and 
the attachment coefficients given by the Gibbs-Thomson 
formula |(5J). Calculations in the left column of Fig. [5] 
are made with 7 = 3.7. The solutions of the Becker- 
Doring equations (black lines) are in a good agreement 
with the results of the kinetic Monte Carlo simulations 



(gray lines), that are repeated from Fig. Q] The coarsen- 
ing laws for adatom and advacancy islands almost coin- 
cide and quickly reach the Lifshitz-Slyozov t 1 / 3 asymp- 
totic. The island size distributions, Fig.^Jc), also almost 
coincide for adatom and advacancy islands, possess ki- 
netic scaling, and agree well with these obtained in the 
kinetic Monte Carlo simulations, cf. Fig.QJd). 

For large bond energies (right column in Fig. [5]), the 
calculations are performed with attachment-limited ki- 
netics, Eq. (fl~5|) , since the kinetic Monte Carlo simula- 
tions point to the Wagner's i 1 / 2 asymptotic. We com- 
pare the discrete distribution of the island energies that 
takes into account the 'magic' island sizes as described in 
Sec. IIIIBl (full black lines) with the continuous one, given 
by the Gibbs-Thomson formula (broken lines). The rela- 
tionship between the discrete and the continuous models 
is established by calculating the energy of a square is- 
land and a circular one with the same number of atoms: 
Eb/k-QT = The calculations are performed for 

7 = 7. The effect of magic sizes is slightly different for 
adatom and advacancy islands. For adatom islands, the 
detachment coefficients b n given by Eq. (jT4]) are excep- 
tionally large for n = m + 1, where m is a magic num- 
ber. Thus, a monomer that has attached to a magic 
island detaches again with a high probability. For adva- 
cancy islands, the detachment coefficients b m for magic 
islands are exceptionally small, so that the detachment 
of an atom from a vacancy island (this atom becomes 
an adatom on the higher level) proceeds at a small rate. 
Both processes make each magic size a trap for further 
island growth, giving rise to the discrete island size distri- 
bution peaked at the magic sizes, shown in Fig. 0] The 
island size distributions presented in Fig. EKc) for this 
case are obtained by averaging over finite ranges of the 
sizes, in the same way as done in the kinetic Monte Carlo 
simulations. 

The time dependence of the average island sizes ob- 
tained for coarsening through the sequence of magic is- 
lands (full black lines in right column of Fig. [5]) are in 
good agreement with the results of kinetic Monte Carlo 
simulations (gray lines). For vacancy islands, the con- 
tinuous island size distribution with the Gibbs-Thomson 
formula gives rise to a notably different coarsening behav- 
ior (broken lines), with a very fast increase of the island 
sizes in an intermediate range. The island size distri- 
butions obtained in the discrete (with magic sizes) and 
the continuous models are also notably different, see Fig. 
03(c). The distribution obtained in the discrete model is 
symmetric with respect to the maximum, similar to the 
one obtained in the Monte Carlo simulations, but notably 
narrower, cf. Fig.QJd). It is worth to note that the distri- 
bution scaled by the average island size does not change 
in time and is the same for the adatom and advacancy 
islands, despite the time evolutions of the average island 
sizes not coinciding and not following a power law. In 
other words, the solution of the Becker-Doring equation 
obeys kinetic scaling in the sense that the island size dis- 
tribution is a function of r/R(t) that does not depend 
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on time. However, R(t) is not described by a power law. 
The continuous model gives a much broader and asym- 
metric island size distribution, shown by broken lines in 
Fig-E^c). The broken-bond counting scheme described in 
Sec. IIII 51 adequately represents the energies E n of small 
islands and quantitatively describes the island size distri- 
bution at the initial stage of coarsening, see Fig. [U How- 
ever, for larger islands it oversimplifies the island energy 
distribution and gives rise to a more narrow distribution 
than found in the simulations. A better model for the 
island energies E n is needed to describe this distribution 
correctly. 

To summarize this section, we show that the Ostwald 
ripening kinetics can be described as an initial value 
problem for ordinary differential equations (|TJ) — dSJ) that 
can be solved by standard numerical methods. This ap- 
proach can be applied to various coarsening problems 
by replacing the Gibbs-Thomson formula ([5]) with Eqs. 
(H3J), ([25]) that admit any dependence of the island en- 
ergy E n on the number of atoms n in it. The alterna- 
tive approach, a numerical implementation of the integro- 
differential equations ([^- (fTTj) ) 68 ! 69 seems much more dif- 
ficult. 



IV. CONCLUSIONS 

Our kinetic Monte Carlo simulations show that the 
Ostwald ripening of 2D islands qualitatively changes 
with increasing bond energy (or decreasing temperature). 



The islands become faceted and the coarsening proceeds 
through a sequence of magic sizes. The Gibbs-Thomson 
chemical potential is not applicable and the detachment 
of monomers from islands is governed by the discrete en- 
ergies of the islands. The coarsening is diffusion limited 
at small bond energies and becomes attachment limited 
at large bond energies. In this latter case, Wagner's 
asymptotic law is reached only after a very long tran- 
sient regime. 

We show that the Becker-Doring equations of nucle- 
ation kinetics are well suited to study the process of 
Ostwald ripening. Two- and three-dimensional coarsen- 
ing processes with diverse limiting mechanisms can be 
simulated by solving a system of ordinary differential 
equations. Concentrations of clusters of all sizes, from 
monomers to ones consisting of millions of atoms, can be 
traced simultaneously. The calculations are not neces- 
sarily based on the Gibbs-Thomson formula but adopt 
any continuous or singular dependence of the cluster en- 
ergy on the number of atoms in it. This approach can be 
applied to a wide range of coarsening problems for two- 
and three-dimensional islands on a surface. 
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